Dynamics of black holes in de Sitter spacetimes 
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Nonlinear dynamics in cosmological backgrounds has the potential to teach us immensely about 
our universe, and also to serve as prototype for nonlinear processes in generic curved spacetimes. 
Here we report on dynamical evolutions of black holes in asymptotically de Sitter spacetimes. We 
focus on the head-on collision of equal mass binaries and for the first time compare analytical 
and perturbative methods with full blown nonlinear simulations. Our results include an accurate 
determination of the merger/scatter transition (consequence of an expanding background) for small 
mass binaries and a test of the Cosmic Censorship conjecture, for large mass binaries. We observe 
that, even starting from small separations, black holes in large mass binaries eventually lose causal 
contact, in agreement with the conjecture. 
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PACS numbers: 04.25.D-, 98.80. Jk,04.20.Dw 



I. INTRODUCTION 

de Sitter spacetimc is the paradigmatic and, in many 
ways, the simplest accelerating universe. It is a maxi- 
mally symmetric solution of Einstein's equations with a 
positive cosmological constant, describing a Fricdmann- 
Robertson- Walker (FRW) cosmology with a constant 
Hubble parameter. Moreover, the large scale structure 
of our universe appears to be that of a de Sitter geom- 
etry, since the large body of observational evidence for 
a present cosmological acceleration is well modeled by a 
positive cosmological constant A [l[ . 

Key questions concerning the evolution towards a de 
Sitter, spatially homogeneous universe are how inhomo- 
geneities develop in time and, in particular ,if they are 
washed away by the cosmological expansion [2J . Answer- 
ing them, requires controlling the imprint of the gravita- 
tional interaction between localized objects on the large- 
scale expansion. Conversely, the cosmological dynamics 
should leave imprints in strong gravitational phenomena 
like primordial black hole (BH) formation 0] or the grav- 
itational radiation emitted in a BH binary coalescence, 
which carry signatures of the cosmological acceleration as 
it travels across the universe. Identifying these signatures 
is not only of conceptual interest but also phenomcnolog- 
ically relevant, in view of the ongoing efforts to directly 
detect gravitational radiation. 

Finally, dynamics in asymptotically de Sitter space- 
times could also teach us about more fundamental ques- 
tions such as Cosmic Censorship: two BHs of sufficiently 



large mass in de Sitter spacetime would, upon merger, 
give rise to too large a BH to fit in its cosmological hori- 
zon. In this case the end state would be a naked singu- 
larity. This possibility begs for a time evolution of such 
a configuration. Does the time evolution of non-singular 
data containing two BHs result in a naked singularity, or 
are potentially offending BHs simply driven away from 
each other by the cosmological expansion? 

Tackling all these issues starts with the study of the 
coalescence process in BH binaries. Since 2005 [||-[(| , 
extraordinary progress has been achieved in the under- 
standing of BH mergers and the associated gravitational 
radiation emission, using numerical methods, while the 
last few years have witnessed a generalized interest in the 
program of extending numerical evolution of BH space- 
times to generic backgrounds [7H15|. In this paper we 
shall take the first step to bring these techniques to a new 
frontier: we report the first numerical evolution of BH 
binaries in an asymptotically de Sitter geometry. Even 
though we consider a range of values for the cosmological 
constant far larger than those which are phcnomcnolog- 
ically viable, these results provide useful insight on the 
general features of dynamical BH processes in spacetimes 
with a cosmological constant, which can improve our un- 
derstanding of our universe. 



II. SCHWARZSCHILD-DE SITTER 

The Schwarzschild-de Sitter spacetime, written in 
static coordinates reads: 
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ds 2 



-f(R)dT 2 + f(R)~ 1 dR 2 + R 2 dn 2 . (1) 
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The solution is characterized by two parameters; the BH 
mass m and the Hubble parameter H and 



f{R) = 1 - 2m/ R - H 2 R 2 , H = ^/IJi 
f(R) has two zeros, at R = R±, i?_ < R + , if 



(2) 



< mH < mH CTit , mH clit = \/l/27 . (3) 

These zeros are the location of the BH event horizon 
(R-) and of a cosmological horizon (R+). If H = 0, then 
R- = 2m; if m = 0, then R + = l/H. If H,m ^ 0, then 
i?_ > 2m and i? + < l/H. Since i? is the areal radius, 
the area of the spatial sections of the cosmological hori- 
zon decreases in the presence of a BH; and the area of the 
spatial sections of the BH horizon increases in the pres- 
ence of a cosmological constant, as one would intuitively 
anticipate. 

The basic dynamics in this spacetime may be inferred 
by looking at radial time-like geodesies. They obey the 
equations (dR/dr) 2 = E 2 — f(R), where r is the proper 
time and E is the conserved quantity associated to the 
Killing vector field d/dT. In the static patch < R < 
R+), E can be regarded as energy. From this equation 
we see that f(R) is an effective potential. This potential 
has a maximum at 
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(m/tf 2 ) 1/3 



(4) 



Geodesies starting from rest (i.e. dR/dr(T = tq) = 0) 
will fall into the BH if i?_ < R < i? max or move away 
from the BH if i? max < R < R + . 

As we will discuss in the next section, the initial data 
for an evolution in the dc Sitter universe can be com- 
puted in a similar manner as has been done in asymptot- 
ically flat space as long as one chooses a foliation with 
extrinsic curvature having only a trace part. Such 
a coordinate system is known for Schwarzschild-dc Sit- 
ter: McVittie coordinates [l6|. These are obtained from 
static coordinates by the transformation (T, R) —> (t, r) 
given by 

* = (i +£ «*K r-M/ /w ™ WJt ,(.) 

where a(t) = exp(Ht) and £ = 2 a(t)r • ^ ne obtains 
McVittic's form for Schwarzschild-de Sitter: 

\] dt 2 + a{t) 2 {l + ^f{dr 2 +r 2 dVL 2 ) . (6) 



For t = constant, one can show that indeed Kj = —H5j. 

By setting m , = in McVittie coordinates one recov- 
ers a FRW cosmological model with k = (flat spa- 
tial curvature) and an exponentially growing scale factor. 
The cosmological horizon He discussed above, located at 
R = l/H, stands at ru c = l/(He Ht ). The spatial sec- 
tions of He seem to be shrinking down in this coordinate 
system. What happens, in fact, is that the exponentially 
fast expansion is taking any observer to the outside of 
He- This is a well known phenomenon in studies of in- 
flation and, as we shall see, has important consequences 
for the numerical evolution. 



III. NUMERICAL SETUP 

The cosmological constant introduces a new term in 
the Hamiltonian constraint obtained after the canonical 
3+1 decomposition: 3 i? - K i} K ij + K 2 = 2A, 
where 3 R denotes the Ricci scalar associated with the 
spatial three metric and K^, K the extrinsic curvature 
and its trace; cf. - In Refs. 0, it was observed 
that imposing a spacetime slicing obeying Kj = —Hdj, 
and a spatial metric of the form dl 2 = ip i jijdx i dx : ' , the 
equations to be solved in order to obtain initial data are 
equivalent to those in vacuum. In particular, for a system 
of N BHs momentarily at rest (with respect to the given 
spatial coordinate patch), the conformal factor ip takes 
the form 



N 



2|r 



(7) 



There are N + 1 asymptotically de Sitter regions, as 
\r — r*(j)| — > 0, +oo; the total mass for observers in 
the common asymptotic region (|r — 77j)| — > +oo) is 

For the numerical implementation we make use of the 
generalized Baumgarte, Shapiro, Shibata and Nakamura 
(BSSN) formulation, e.g. [17|,l20|. For the case of van- 
ishing cosmological constant the definition of the BSSN 
variables %, 7y, Aij, T k and their evolution in time is 
given by Eqs. (1, A1-A8) in Ref. [MQ. For our simu- 
lations with A / we apply two modifications to this 
formalism, (i) The evolution equation for K becomes 
(d t — Cp) K = [•••] — aA, where [. . . 1 de notes the "stan- 
dard" right-hand side of Eq. (A6) in [20|. (ii) a new vari- 
able x — cxp(2Ht)x has been evolved instead of x 0- 
The reason is that for BH evolutions it is crucial to im- 
pose a floor value on x, typically 10~ 4 or 10~ 6 , which is 
inconsistent with the natural behaviour of this variable 
in a de Sitter spacetime: x _1 ~ exp(2iJ<). In contrast 
X — > 1 when r — > oo for all times. The evolution equation 
for x is given by 



d t x = 2x(aK - d^/'d + (3 l d t x + 2H X , 



(8) 



and replaces Eq. (A6) of [20|. Boundary conditions for 
all quantities are imposed by looking at the behavior of 
massless perturbing fields in a pure de Sitter background. 
Accordingly, we impose the following asymptotic behav- 
ior for all BSSN variables 

dtf-d t h + -±-d r f+ f —^-H(f-f a )=Q. (9) 
a(t) a(t)r 

We should note that we also performed evolutions using 
different sets of boundary conditions, to test the indepen- 
dence of the results on boundary conditions imposed in a 



1 We note a missing factor of \ m the final term of Eq. (A6) therein. 
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region with no causal contact with the interaction region. 
As far as the behaviour and location of the horizons and 
all quantities discussed in this paper are concerned, no 
noticeable difference could be found. 

Our numerical simulations use the Lean code 20 1 
which is based on the Cactus computational toolkit 21 1 
and the Carpet mesh refinement package [22], [23| . The 
calculation of Black hole Apparent Horizons (BAHs) and 
Cosmological Apparent Horizons (CAHs) is performed 
with AHFinderDirect [U,[25|. We remark that BAHs, 
found as marginally trapped surfaces, indicate in the Sit- 
ter space (with the same legitimacy as in asymptotically 
fiat space) the existence of an event horizon |26| ■ CAHs 
are surfaces of zero expansion for ingoing null geodesies. 
In a single BH case, in McVittie coordinates, the BH 
event horizon and cosmological horizon are indeed foli- 
ated by apparent horizons. 

The "expanding" behaviour of the coordinate system 
led us to add a new innermost refinement level at periodic 
time intervals so as to keep the number of points inside 
the cosmological horizon approximately unchanged. The 
necessity for adding extra refinement levels effectively 
limits our ability to follow the evolution on very long 
timescales, as the number of time steps to cover a fixed 
portion of physical time grows exponentially. This fea- 
ture resembles in many ways the recently reported work 
by Pretorius and Lehner on the follow-up of the black 
string instability p^ . 



IV. NUMERICAL RESULTS 




FIG. 1. Conformal factor \ for a single BH evolution with 
H = 0.8-Hcrit using the McVittie slicing condition, Eq. (|10[1 , 
The obtained numerical results are plotted, along the z coor- 
dinate (symmetry \{— z) = x( z ) imposed at z = 0), against 
the expected analytical solutions (solid lines). 



As a first test on the numerical implementation, we 
performed evolutions of a single BH imposing the McVit- 
tie slicing condition; that is, we use ([5]) as initial data and 



impose 

d t a = 4mrHe m /(m + 2re m ) 2 , d t f3 l =0, (10) 

throughout the evolution. The analytical solution © can 
be compared with the numerical results. For a single BH 
evolution with m = 1 and H = 0.8iJ C rit, the results are 
displayed in Fig.[T] Using this slicing, the runs eventually 
crash (at t ~ 12m). By contrast, the standard "1+log" 
slicing condition 



d t a = ftdiU -2a(K - K ) 



[11) 



where K = —3H = — v3A, enables us to have long 
term stable evolutions. As consistency checks, the areal 
radii at the apparent horizons (both BH horizon and cos- 
mological horizon) are constants in time and have the 
value expected from the analytical solution in a single 
BH spacetime. Moreover, the areal radius at fixed coor- 
dinate radius evolves with time in the way expected from 
the exact solution. 

For binary BH initial data, we start by reproducing 
the results of Nakao et. al [19| . where the critical dis- 
tance between two BHs for the existence of a common 
BAH already at t = was studied. We thus prepare 
initial data ([7]) with iiij = m 2 & n d take all quantities in 
units of the total mass m = mi + m 2 . The two punc- 
tures are set initially at symmetric positions along the 
z axis. The critical value for the cosmological constant, 
for which the BH and cosmological horizon coincide is 
now mH CT it = l/y/27. We call small (large) mass bi- 
naries those, for which H < H cri t (H > -f/crit)- Our 
results for the critical separation in small mass binaries, 
at t = 0, as function of the Hubble parameter are shown 
in Fig. [3J The line (diamond symbols) agrees, after a 
necessary normalization, with Fig. 14 of [19j . 

We now consider head-on collisions of two BHs with 
no initial momentum, i.e. the time evolution of these 
data. We have monitored the Hamiltonian constraint 
violation level for cases with and without cosmological 
constant. We observe that the constraint violations are 
comparable in the two cases and plot in Fig. [3] a snap- 
shot of the Hamiltonian constraint violation at t = 48m 
for parameters H = 0.9if C rit and d = 0.8m, a typical 
case with non-zero cosmological constant. We have used 
two resolutions, m/160 and m/192 (on the innermost re- 
finement level) and have rescaled the dashed curve by 
Q2 = (192/ 160) 2 as expected for second-order conver- 
gence. 

For subcritical Hubble constant H < -ff C rit = 
l/(\/27m), we monitor the evolution of the areal ra- 
dius of the BAHs and that of the CAH of an observer 
at z = 0. For instance, for H = 0.9H CI a and proper (ini- 
tial) separation 3.69m we find that the areal radii of the 
BAH and CAH are approximately constant and equal to 
i?BAH — 2.36m and Rcak — 4.16m, respectively. As ex- 
pected the two initial BAHs, as well as the final horizon, 
are inside the CAH. As a comparison, a Schwarzschild-dc 
Sitter spacetime with the same H has -Rbah — 2.43m and 
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• • BAH forms 
x — x no BAH forms 
- geodesic prediction, m =0.7 
i=0 




FIG. 2. Critical coordinate distance for small mass binaries, 
from both initial data and dynamical evolutions, as well as 
a point particle estimate, as a function of H/H cl i t . We ob- 
tain this estimate from the coordinate distance to the horizon, 
Eq. (|4]), for a particular value of m. The t = line refers to 
the critical separation between having or not having a com- 
mon BAH in the initial data. The inset shows details of the 
approach to the critical line for H = 0.6i? cr it, where a is an 
acceleration parameter. 
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FIG. 3. Hamiltonian constraint violation along the z-axis at 
time t = 48m for a simulation with H = 0.9_ff C rit and initial 
distance d = 0.8m. 



i?CAH — 4.16771. This suggests that the interaction effects 
(binding energy and emission of gravitational radiation) 
are of the order of a few per cent for this configuration. 

As the initial separation grows, so does the total time 
for merger. For separations larger than a critical value, 
the two BHs do not merge, but scatter to infinity. For 
such scattering configurations, the simulations eventually 
exhibit a regime of exponentially increasing proper dis- 
tance between the BAH. Just as in scatters of high energy 
BHs 28[, here we find that the immediate merger/scatter 
regimes are separated by a blurred region, where the 
holes sit at an almost fixed proper distance for some 
time; cf. Fig. 2J By performing a large set of simula- 



tions for various cosmological parameters H and initial 
distance d, we have bracketed the critical distance for 
the merger/scatter region as a function of the Hubble 
parameter H for the "dynamical" case, i.e., the initial 
coordinate distance between the BHs such that no com- 
mon BAH forms. The results are displayed in Fig. [5] 
(circles and x symbols). 
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FIG. 4. Proper distance between the BH horizons as a func- 
tion of time for the H = 0.9// C rit, and initial (coordinate) 
distance d ~ 0.9m. The two holes stay at approximately 
constant distance up to t « 8m after which cosmological ex- 
pansion starts dominating. 

As expected the critical distance becomes larger as 
compared to the initial data value ("i = 0" line): there 
are configurations for which a common BAH is absent in 
the initial data but appears during the evolution (just as 
in asymptotically flat spacetime). The numerical results 
can be qualitatively well approximated by a point par- 
ticle prediction - from Eq. [4] To do such comparison a 
transformation to McVittie coordinates needs to be done; 
we have performed such transformation at McVittie time 
t = 0. Intriguingly, for a particular value of m ~ 0.7, 
the point particle approximation matches quantitatively 
very well the numerical result; the curve obtained from 
the geodesic prediction in Fig. [3] is barely distinguishable 
from the numerical results. 

A further interesting feature concerns the approach to 
the critical line. For an initially static binary close to the 
critical initial separation, the coordinate distance d scales 
as d = do + at 2 . In general the acceleration parameter 
scales as log a = C + riog(ci — do), where T = 1 in the 
geodesic approximation. A fit to our numerical results 
for H = 0.6-ff cr it (dashed curve in the inset of Fig. [2]) for 
example yields C = —3.1, T = 0.9 in rough agreement 
with this expectation. Details of this regime are given in 
the inset of Fig. [5] 

Finally, we have performed evolutions with H > H crl f 
On the assumption of weak gravitational wave release, 
such evolutions can test the Cosmic Censorship conjec- 
ture since the observation of a merger in such case would 
reveal a violation of the conjecture [2{|. From general 
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FIG. 5. Snapshots at different times (from left to right t/m = 
0.0,8.0156,20.016) of a simulation with H = 1.05-Hcrit, and 
an initial coordinate distance d/m = 1.5002. The dotted blue 
line denotes the CAH (for an observer at z = 0) which is first 
seen in this simulation at t/m = 8.0156, highly distorted. At 
late times, the CAH has an areal radius of R — 4.94876, while 
the "theoretical value" for pure dS is R = 1/H = 4.94872, a 
remarkable agreement showing spacetime is accurately empty 
dS for the observer at z = and the BHs are not in causal 
contact. 



control. Our results open the door to new studies of 
strong field gravity in cosmologically interesting scenar- 
ios. In closing, we would like to mention that our results 
are compatible with Cosmic Censorship in cosmological 
backgrounds. However, an analytic solution with multi- 
ple (charged and extremal) BHs in asymptotically dc Sit- 
ter spacetime is known, and has been used to study Cos- 
mic Censorship violations [30j . In collapsing universes a 
potential violation of the conjecture has been reported, 
although the conclusion relied on singular initial data. To 
clarify this issue, it would be of great interest to perform 
numerical evolution of large mass BH binaries, analogous 
to those performed herein, but in collapsing universes. 
This will require adaptations of our setup, since the "ex- 
panding" behaviour discussed of the coordinate system 
will turn into a "collapsing" one, which raises new nu- 
merical challenges. 



arguments and from the simulations with H < H CT it> we 
know the cosmological repulsion will dominate for suf- 
ficiently large initial distance and in that case we can 
even expect that a CAH for the observer at z — will 
not encompass the BAHs. This indicates the BHs are no 
longer in causal contact and therefore can never merge. 
Our numerical results confirm this overall picture. To 
test the potentially dangerous configurations, we focus 
on the regime in which the BHs are initially very close. 
A typical example is depicted in Fig. [5j for a supercrit- 
ical cosmological constant H = 1.05-ff cr it, and an initial 
coordinate distance d/m = 1.5002. Even though the ini- 
tial separation is very small, we find that the holes move 
away from each other, with a proper separation increas- 
ing as the simulation progresses. In fact, further into the 
evolution, a distorted CAH appears, and remains for as 
long as the simulation lasts. At late times, this CAH 
is spherically symmetric, and has an areal radius which 
agrees, to within 10 -5 , with that of an empty de Sitter 
spacetime with the same cosmological constant. The evo- 
lution therefore indicates that the spacetime becomes, to 
an excellent approximation, empty dc Sitter space for the 
observer at z = and that the BHs are not in causal con- 
tact. Observe that qualitatively similar evolutions can be 
found in small mass binaries when the initial distance is 
larger than the critical value 

V. FINAL REMARKS 

We have presented evidence that the numerical evo- 
lution of BH spacetimes in de Sitter universes is under 
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